Astronomy & Astrophysics manuscript no. aa'2003'0950 


February 2, 2008 


(DOI: will be inserted by hand later) 





Analytic Solutions for the Evolution 
of Radiative Supernova Remnants 



R. Bandiera 1 and O. Petruk 2 



o 
o 

(N 

in 

(N 



> 
oo 
o\ 
in 

(N 

o 

o 
^ ' 

6 ■ 
a: 



1 Osservatorio Astrofisico di Arcetri, Largo E.Fermi 5 Firenze 50125, Italy 
e-mail: bandiera@arcetri . astro . it 

2 Osservatorio Astrofisico di Arcetri, Largo E.Fermi 5 Firenze 50125, Italy 

Institute for Applied Problems in Mechanics and Mathematics, Naukova St. 3-b Lviv 79000, Ukraine 
e-mail: petruk@arcetri . astro . it 

Received 7 January 2004 / Accepted 10 February 2004 

Abstract. We present the general analytic solution for the evolution of radiative supernova remnants in a uniform interstellar 
medium, under thin-shell approximation. This approximation is shown to be very accurate approach to this task. For a given 
set of parameters, our solution closely matches the results of numerical models, showing a transient in which the deceleration 
parameter reaches a maximum value of 0.33, followed by a slow convergence to the asymptotic value 2/7. Oort (1951) and 
McKee and Ostriker (1977) analytic solutions are discussed, as special cases of the general solution we have found. 
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1. Introduction 

In recent years numerical modelling of the structure and evo- 
lution of supernova remnants (SNRs) has reached an unprece- 
dented level of detail. Nonetheless analytic models still play 
a very important role, when general properties have to be in- 
vestigated, as well as when direct relations have to be drawn 
between pure observational quantities (like size, flux, etc.) and 
intrinsic physical parameters. 

The adiabatic phase of SNR evolution (in a uniform and 
homogeneous medium) is well described by the Sedov (1959 1 
analytic solution, which reproduces both the SNR radial evolu- 
tion and its inner structure. This exact solution has been made 
possible by the fact that during this phase the SNR evolution 
is self-similar. This is no longer the case when radiative losses 
become important, and therefore no exact analytic solution is 
known for the late SNR evolution. 

Approximated solutions in the adiabatic regime and be- 
yond may be also obtained using a "thin-shell" model (see e.g. 
Zel'dovich & Raizer lT966l Ostriker & McKee IT9881 . This ap- 
proach assumes that the whole mass (and therefore kinetic en- 
ergy) of the SNR is located in a rather thin shell just behind 
the outer shock; while the inner region is filled with a very hot 
and rarefied gas, of negligible total mass, but containing most 
of the SNR internal energy. 

For the adiabatic phase this approximation is only moder- 
ately accurate (see e.g. Zel'dovich & Raizer 1966 1. In fact, ac- 
cording to the Sedov solution the gas density vanishes in the 



inner regions while its pressure keeps finite; however, the outer 
layer containing most of the mass is geometrically rather thick. 
On the other hand, numerical works (e.g. Falle 1975 Blondin et 
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al.[T998, hereafter BWBR) trace the formation of a much thin- 
ner shell in the radiative phase, therefore indicating that a thin- 
shell approximation should be far more accurate in describing 
the late evolution. 

Oort dl95U presented a first thin-shell approach to a radia- 
tive SNR expansion. By assuming momentum conservation in 
the shell, he found the SNR radius to evolve as R cc t 1 ^. This 
solution, also known as "momentum-conserving snowplow", 
assumes that cooling is extremely efficient everywhere (and 
therefore that the interior pressure vanishes). However, numer- 
ical models (e.g. Chevalier 1974 1 show that, even in the radia- 
tive phase, the gas in the central regions becomes so rarefied 
that its cooling time still keeps considerably longer than the 
SNR age. This led McKee & Ostriker J1977I to introduce a 
"pressure-driven snowplow" model, in which a fossil pressure 
in the hot interior has a substantial dynamical effect on the outer 
shell: in this case the radial evolution is R cc t 2 ^ (for adiabatic 
index y = 5/3). 

Even though the "pressure-driven snowplow" formula gets 
closer than the "momentum-conserving snowplow" one to the 
numerical results, some discrepancy still remains. For instance, 
by defining the "deceleration parameter" as m = d log R jd log f, 
numerical models obtain an asymptotic value ranging from 
0.31 (Chevalier [19741 to 0.33 (BWBR). These values are sig- 
nificantly different from the analytic value, 2/7 (namely 0.286), 
and various authors have discussed the origin of such discrep- 
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ancy. Cioffi et al. ( 1988 1 ascribe it to a "memory" of the previ- 
ous Sedov phase, leading to an actual internal pressure larger 
than that derived from the analytic model. BWBR, instead, at- 
tribute this discrepancy to the influence of the reverse shock, 
which moves towards the center raising the thermal energy, 
thus leading to a milder deceleration. 

Other authors have estimated analytically the radial evolu- 
tion under more general conditions than those given above. 

Ostriker & McKee ( 1988 1 have shown that, for a general y 
as well as a power-law ambient density profile (p a (r) °c r~ w ), 
m = 1/(4 — cj) for a "momentum-conserving snowplow", while 
m = 2/(2 + 3y - a>) for a "pressure-driven snowplow". 



Liang & Keilty (2000 1 have considered the case in which 
only a (constant) fraction e of the kinetic energy of the in- 
coming flow is radiated in the outer shock. For y = 5/3, m 
is found to decrease quasi-linearly with e, from 2/5 for the adi- 
abatic case (e = 0) to 2/7 for the fully radiative case (e = 1); 
and a value of e of about 0.8 (0.6) is required in order to ob- 
tain m = 0.31 (0.33), as indicated by the numerical models. 
However, while e < 1 may be appropriate to describe gamma- 
ray burst afterglows (Cohen et al. 1998 1, SNR radiative shocks 
should be described as fully radiative ones (namely with e very 
close to unity). 

The effect of cooling in the hot interior on the deceleration 
parameter has been studied by Gaffet ( 1983), with the follow- 
ing results. Adiabaticity holds throughout most of the volume 
occupied by the hot gas, while cooling occurs only near the 
boundary with the radiative shell, giving as effect a net mass 
transfer from the hot interior to the shell. Assuming that the 
gas in the hot interior follows a cooling law A oc T~°, this pa- 
per discusses different regimes for different choices of y and c, 
showing that the asymptotic value of m must be in the range 
between the values 1 /4 (Oort limit) and 2/(2 + 3y) (McKee and 
Ostriker limit). 

When c > (2/3)(5y - 8)/(3 - 2y) (c > -2/3 for 7 = 5/3) 
the asymptotic value of m is 1/4. This is the case for cooling 
functions typical of the SNR regime, where c > and is usu- 
ally taken in the range from 0.5 (e.g. Cioffi et al. 1 19881 1 to 1.0 
(BWBR). Therefore, according to this result, after the com- 
monly known radiative phase there could be a very late evo- 
lutive phase, in which radiative losses of hot interior begin to 
be prominent and the SNR evolves as "momentum-conserving 
snowplow". However, as it may be derived from numerical re- 
sults, in a typical SNR, radiative cooling of the hot interior 
is negligible until very late times. Therefore the onset of the 
"momentum-conserving snowplow" regime should occur only 
near the end of a SNR lifetime, or not to occur at all. 

A common limitation of all above-mentioned analytical 
models is that the radial evolution of the radiative shock has 
been approximated by a power-law behaviour R oc f" (with 
constant m). This allows a simplified treatment of radiative 
SNR evolution; however, it is natural to expect that a power- 
law expansion occurs only at late times (i.e. at large R values), 
after the transition from adiabatic to radiative expansion has 
been completed. 

In this paper we shall show: 1) that the quoted difference 
between the numerical and asymptotic analytic value is just 
a consequence of the fact that the time needed to reach the 



asymptotic power-law regime is long compared with the age 
of the SNR; 2) that the SNR radial evolution during the tran- 
sient phase is adequately described by a thin-shell model; 3) 
that a general analytic solution of this problem exists. 

The plan of the paper is as follows. In Sect. 2 we derive 
a differential equation that describes the radiative SNR evolu- 
tion, find its complete solution for an arbitrary adiabatic coef- 
ficient, discuss the general properties for the two branches of 
solutions that we find, and show that the Oort solution is just a 
special case of our general solution. Sect. 3 focuses on the con- 
ventional case y = 5/3, which allows simpler analytic formu- 
lae, and shows that our solution tends to the asymtotic regime 
given by the solution of McKee & Ostriker d!977> : with the use 
of numerical results, we also derive the most appropriate initial 
conditions for our solution. Sect. 4 concludes. 



2. Equations and solutions for a general adiabatic 
index 

Let us consider a fully radiative shock expanding into a uni- 
form medium and neglect the cooling of the hot interior. In the 
thin-shell approximation, mass, momentum, and central pres- 
sure evolution are described by the following set of equations: 



dM 
dt 



= 4np a R 2 R, 



— = 4nPR 2 , 

dt 

dP „ R 

— = ~3yP~, 
dt ' R 



(1) 
(2) 
(3) 



where R and R are respectively shock radius and velocity, M is 
the mass of the radiative shell, P is the pressure of the (adiabat- 
ically evolving) inner region, y is the adiabatic coefficient, and 
p a is the (constant) density of the ambient medium. 
The above equations can be reduced to a single one: 



R+ — = 3P ° R ° 7 R -3v-i 
R Pa 



(4) 



where the quantities R and P indicate respectively the SNR 
radius and the pressure of the hot cavity at a reference time (? D ), 
that can be arbitrarily chosen. 

In order to solve analytically Eq. (0J for a general value of 
y (with the condition 1 < y < 2) let us first define the quantity: 



K 



2 PR?y 

2-7 Pa 



(5) 



K is constant in time, and therefore it can be evaluated in terms 
of quantities at the time t Q . By using the substitution w(R) — R 2 , 
Eq. (|4) trasforms into: 



(6) 



that is a linear differential equation and can then be easily inte- 
grated. Its general solution is: 



K(R 



-3y . 



HR\ 



(7) 
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where the constant H is, at any time, equal to: 
// = /?3(2-r)( 1 _ 



(2 - y)pa , . 2 



2P 



(8) 



In particular, it may be expressed in terms of quantities at the 
time t Q . Depending on the sign of H, there are two different 
branches of solutions. By evaluating the kinetic energy of the 
shell and the thermal energy of the inner hot bubble respec- 
tively as: 



^kin — 



Eth = 



4nR i p. d R z 



3 

4nR 3 



2 
P 



3 y- 1' 



(9) 



(10) 



Eq. l|8} shows that, in the two branches, the energy ratio x = 
Ekm/Eth = (y-l)p a ^ 2 /2-P is respectively less (//-positive case) 
and greater (//-negative case) than (y - 1)/(2 - y), and that time 
evolution does not change the sign of this inequality. 

Let us label these two branches of solutions as "slow" 
and "fast", depending whether the kinetic energy is respec- 
tively less (H positive case) or greater (H negative case) than 
(y- l)/Sth/(2-y); or, equivalently, less or greater than (y-l)E tot 
(where E tot = E^ n + E±). The choice of the appropriate branch 
of solutions only depends on the initial conditions. 

Although in the next section we shall see that the slow case 
is that physically relevant for the SNR evolution, let us discuss 
here both branches. When double signs are shown, in some of 
the following equations, the convention used is that the upper 
sign refers to the slow branch, while the lower sign to the fast 
branch. Once defined space and time scale units as 

R = (±//) I/3(2 - y) ; (11) 
t = (±/f ) (2+3r)/6(2-r) A: -i/2 i (12) 

and introduced the dimensionless space and time coordinates 
r = R/R, t = t/t, the evolution in size follows the equation: 



— = Vr- 3 r + r- 6 , (13) 
dr 

that can be integrated to give r(r). The dimensional velocity is 
obtained multiplying dr/dr by the velocity scale V = R/t. It is 
evident that, while fast solutions extend to all positive values of 
r, solutions in the slow branch are real only for r > 1 . 

For a general value of y, the solution involves hypergeo- 
metric functions (F), and can be written as: 



r F (r) = ^-F\\, 



,1 + 



4 \T 3(2-y)' 3(2- r )' 



■ _ r 3(2-y) 



1 + 



4 \2' 3(2-7)' 3(2-7)' 
for the fast branch, and as: 



l\ + C 



(14) 



ir 4 \ 4 

r s (r) = — F -, . . 

sw 4 \2' 3(2-7) 3(2-7) 



l + ^_ ;r 3(2- y) 



-MI. 



i + 



1 +C 



(15) 



4 \2' 3(2-7)' 3(2-7)' 

for the slow branch. The time evolution of the SNR radius is 
obtained by inverting the above equations. Note that individual 



terms in Eq. (1151 are complex, but when r > 1 their imaginary 
parts cancel out. Eqs. dl4> and dl5> contain an arbitrary con- 
stant, C; both equations have been written here in such a way 
thatC = r(l). 

A quantity useful to describe the evolution is the decelera- 
tion parameter m. The general formula for this quantity is rather 
complex, but its asymptotic behaviour at large values of r may 
be evaluated as: 

2 , 6(2-7) 

m(r) — 



_ r -3(2-r) + 



0( r -6(2-7)).(16) 



2 + 37 (97 - 10)(2 + 37) 
This power expansion is valid in the range 10/9 < 7 < 2. 
Note that the limit 2/(2 + 37) is the same found by Ostriker and 
McKee fTWgl . 

In terms of the dimensionless quantity r, the ratio of kinetic 
and thermal energies is: 



2 - 7 v 



-3(2-7)1 



(17) 



Therefore at late times the asymptotic value of this ratio is x as = 
(7 - l)/(2 - 7) for both branches. 

The ratio between internal and shock pressure is: 

P_ = (7-IX7+D = (2-7X7+I) 
P, 4x 4(1 +r- 3 ( 2 -r))' 

then leading to the asymptotic value (2 - 7X7 + l)/4 
The total energy follows the evolutive law: 

Ep 



(18) 



, (19) 
\7~ 1 

where E = InpJ/ 2 ^ /3. It is easy to show that, in the range 
of validity of the solutions, dE tot /dr is always negative, as ex- 
pected for a radiative solution. 

Finally, it can be shown that the original Oort J195U solu- 
tion, R oc f 1 / 4 , is just a special case of the fast-branch solution. 
Neglecting the pressure effect leads the right side of Eq. ® to 
vanish. Since P = implies the energy ratio x to diverge (and 
therefore to be larger than (7 - l)/(2 - 7)), the solution must 
belong to the fast branch. From Eq. it is apparent that K 
is required to vanish, while H — > -00, in such a way that the 
product -KH be equal to R 6 R 2 (being this a constant, it can be 
then evaluated in terms of R and /? ). Using Eqs. dl It and (I12> . 
it can be shown that both R and I diverge, so that this solution 
must be limited to vanishingly small r and r values. Therefore 
Eq. Jl 31 simplifies into dr/dr - r~ 3 , which admits the solution 
r(r) = r 4 /4-l/4 + C, where we have defined C = t(1) for con- 
sistency with the formulation given in Eq. dl4> . For C — 1/4, 
which means r(0) = 0, we simply have r(r) = r 4 /4 that, when 
inverted, gives the Oort's law r oc r 1 ^ 4 . The same result can be 
extracted from the general solution, Eq. J14l> . by using the fact 
that F(a, b, c; x) — > 1 when x — > 0. 

3. The slow branch of solutions for 7 = 5/3 

In the standard case y = 5/3, Eqs. dl4ll and dl5l l get a much 
simpler functional dependence, respectively: 

tfW 



? 1 Q a/o 

— V7+T(5r 3 - 6r 2 + 8r- 16) + — + C, 

35 35 



Ts (r) = — Vr- l(5r 3 + 6r 2 + 8r + 16) + C, 



(20) 
(21) 
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where again we use C = t(1). 

The deceleration parameter m, for the two branches, evalu- 
ates: 



0.4 



mp(r) 



m s (r) = 



2 r+l 
35~ 



-(5r-6r 2 + Sr- 16) + 



V^ + T/ 18 



r 4 



( c+ §4 (22) 



2 r- 1 
35 "T 4 " 



-(5^ + 6r z + 8r + 16) + , C. (23) 



The asymptotic behaviour at large values of r is: 



2 2 / 1 



(24) 



(which is consistent with the more general Eq. dl6» . Therefore, 
for r approaching to infinity, in both branches m tends to the 
value 2/7, namely to the asymptotic solution given by McKee 
and Ostriker I1977> . However, an analysis of Eqs. i22i and 
j23l > show different properties for the two branches. In partic- 
ular, only in the slow branch m(r) shows a local maximum. In 
a given solution, the position of the maximum and the value 
reached by m are related by: 



2(r„ 



1) 



(25) 



valid for r max > 8/7. Therefore, ;« max is always larger than 2/7 
for any solution in the slow branch, and it can be considerably 
larger than 2/7, if r max gets close to 8/7. Furthermore, if during 
its evolution m(r) is larger than 2/7 and still increasing with r, 
it must reach a maximum before approaching the asymptotic 
value 2/7, and then it must belong to the slow branch. This 
is what shown by BWBR (their Fig. 3, actually limited to the 
increasing part): therefore in the following we shall consider 
only the slow branch of solutions. 

Let us use the numerical results from BWBR to deter- 
mine the most appropriate parameters for our analytic solution. 
For the fit we use the evolution of m (BWBR, their Fig. 3), 
excluding the oscillatory transient: the numerical data fitted 
are for times ranging from 7.4 x 10 4 till 3.0 x 10 5 yr. The 
best (least-square) fit is obtained for C = -0.248 ± 0.006; 
while the time scale inferred from this fit allows us to fix 
t = (3.64 ± 0.05) x 10 4 yr (for all best fit quantities, here we also 
indicate their l-cr error). In Fig.^ the best fit curve is shown 
against the numerical data. The best fit curve reaches its maxi- 
mum value (ra max = 0.328) at r max = 2.11 (i.e. at r max = 6.18). 
Moreover, using Fig. 8 (velocity evolution) from the BWBR, 
we derive R = 17.6 + 0.1 pc. From these quantities, the dimen- 
sional scaling for energy is E = (1.51 ± 0.07) x 10 51 erg. 

The numerical simulation we refer to corresponds to the 
following basic physical parameters: energy of the explosion, 
£gi = 1 (in units of 10 51 erg); and hydrogen ambient density, 
n a = 0.84 cm 3 . Adopting the same definitions of BWBR for 
the transition time, and corresponding SNR radius 



t trdn * 2.9xl0 4 £ 4 ( 1 V /17 



yr. 



19.1<V /17 PC 



(26) 
(27) 




i ill i i 



4.6 4.8 5 5.2 

log 10 (t/yr) 

Fig. 1. Our best fit solution (slow branch, C = -0.248, t = 
3.64 x 10 4 yr) compared with numerical data from BWBR. The 
small frame shows just the data used for the fit. With the ex- 
ception of strong oscillations in the early transient, the analytic 
solution closely describes also the evolution at earlier times. 
Dashed lines indicate the positions of f max and fi nters , while the 
dotted line is obtained by using the analytic fit by Cioffi et al. 
1988 (see text). 



we determine the dimensionless quantities ?/f tran = 1.14 and 
R/Rtmn - 0.85. All these quantities, although obtained after the 
comparison with a specific numerical simulation, can be taken 
of general validity, for a SNR expanding in a homogeneous 
medium, because the analytic solutions allow scaling. 

Fig. 1 shows also (dotted curve) the analytic fit as derived 
from Cioffi et al. fHSggl l. Cioffi et al. fi^SSl use a different cool- 
ing function from that of BWBR, whose simulation is shown 



in Fig. 1. Therefore in Cioffi et al. ([T988) the functional de- 
pendence of time and length scales on the model parameters 
are different from those given in Eqs. \26\ and j27l >. In order 
to compare our results with Cioffi et al. (Q588 1 analytic fit, we 
have evaluated the scaling time f tran (labelled as fpos in their 
Eqs. (3. 10) and (3. 11)), using E 51 = 1 andn a = 0.84 cm" 3 (with 
solar abundances), obtaining f tran = 14670 yr. The formula used 
for the deceleration parameter is m = 0.3/(1 - t lrdn /4t), as de- 
rived from Eqs. (3.32) and (3.33) in that paper. It is apparent, 
from Fig. 1, that Cioffi et al. ( 1988) fit does not trace the evo- 
lution of the deceleration parameter m. 

Since the evolution of the SNR radius is a continuous func- 
tion of time, let us compute the time at which the radiative 
solution intersects the Sedov one. With the parameters given 
above, it happens at f inters = 1.16f tran (when fl inters = 1.06R aim ), 
namely, using our dimensionless variables, at Ti nte rs = 1-01 
(with r inters = 1.24). At this time, the SNR kinetic, Eq. 0, 
and thermal, Eq. fllOi . energies are respectively E^ n = 0.191 x 
10 51 erg and £th = 0.489 x 10 51 erg, equivalent to a total energy 
E tot - 0.680 x 10 51 erg, and to an energy ratio x - 0.390. 

Nicely, although fortuitously, at fj ntels the value of t is very 
close to unity, while that of x is very close to the Sedov value 
(^Sed = 0.394). We could then use t q = 1 and x = ^s e d as an 
approximate criterion, from which to derive, analytically, that: 
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1) the solution must belong to the slow branch, since x a < % as ; 

2) r Q and m Q , evaluated using the relatioships: 

r = 2/(2 - x ), (28) 
m Q = r (2-^ ) 7/ X /2 /16, (29) 

are respectively r Q a 1.245 and m a 0.206; 3) C, evaluated 
using Eq. J21l i. is a -0.272 (to be compared with the best fit 
value -0.248). 

4. Conclusions 

Using a thin-shell approach, we have developed the definitive 
analytic treatment for the evolution of a SNR in the radiative 
phase, and we have also obtained a series of interesting rela- 
tions. The main findings of the present work are the following. 

The discrepancy between the analytic prediction of the 
asymptotic value of the deceleration parameter (m = 2/7, 
McKee and Ostriker ri977> and that derived numerically (m = 
0.33, BWBR) is only apparent. This discrepancy has been at- 
tributed to the presence of a reverse shock moving towards the 
center. We show, instead, that a thin-shell model, that by defini- 
tion does not contain any information on inner structure details, 
closely fits the SNR evolution as derived numerically. 

We then confirm that 2/7 is the correct asymptotic value, 
even though the convergence towards this value is expected to 
be slow. We believe that, if BWBR numerical simulation had 
been runned until later stages of the SNR evolution, it would 
have shown that m does not keep constant to 0.33, but eventu- 
ally decreases. This has been already pointed out by Chevalier 
SW]4\ and can be seen in Fig. 5 of Cioffi et al. ( fT988b . in Fig. 3 
of Falle (HUB, and in Fig. 4b of Mansfield & Salpeter JT9741 . 
However, the convergence to the asymptotic value may need 
times longer than the SNR lifetime. 

It might be expected that the evolution will eventually 
change from a "pressure-driven snowplow" (m = 2/7, McKee 
and Ostriker 119771 to a "momentum-conserving snowplow" 
(m = 1/4, Oort 1951), as a consequence that the right side 
of Eq. @ vanishes when R — > oo. However, Cioffi et al. ( 1988 ) 
have noted that, even at very late times (~ 10 2 f tran ), the decel- 
eration parameter m is still closer to 2/7 than to 1/4. We have 
shown that such evolutive transition may in fact not occur, be- 
cause the two kinds of evolution are associated with two dif- 
ferent branches of solutions, corresponding to different initial 
conditions. In other words, for a deceleration parameter smaller 
than 2/7, the right side of Eq. @ (with y = 5/3) vanishes more 
slowly than the left side. Therefore, unless the pressure term 
is negligible from the beginning, or part of the internal energy 
of the hot interior is lost by other channels (e.g. by electron 
conduction, or radiative processes), pressure effects must play 
an important role in the evolution at any time, until the SNR 
merges with the ambient medium. 

Some of the conclusions we have presented here could have 
been reached long before. In fact, Blinnikov et al. (|1982 ) have 
reduced the system of equations for the thin-shell evolution to 
the single equation, equivalent to our Eq. @, and have obtained 
a solution equivalent to our Eq. Mil . However, they have not 
discussed the properties of this solution, taking that it would 
have quickly relaxed to the asymptotic behaviour. 
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